Alpha diversity

Paths and libraries setting

# Load main packages, paths and custom functions
source("../../../source/main_packages.R")
source("../../../source/functions.R")

# Load supplementary packages
packages <- c("multcomp", "openxlsx", "kableExtra", "dplyr")
invisible(lapply(packages, require, character.only = TRUE))

Load phyloseq object after decontam

ps.filter <- readRDS("../../../../output/1_MED/1D/1D_MED_phyloseq_decontam.rds")

Alpha diversity plots

Whole

# select whole samples
ps.whole <- subset_samples(ps.filter, Organ=="Whole")
ps.whole@sam_data %>% data.frame() -> test

# estimate alpha diversity

measures <- c("Observed", "Chao1", "Shannon")

p1 <- plot_richness(ps.whole, 
                      x="Sample", 
                      color="Strain", 
                      measures=measures, 
                      nrow = 1)

# extract data from richness plot to custom it
df <- p1$data

# change levels and order of species and location for plot
levels(df$Species) <- c("AA", "CP", "CQ")
df$Species <- factor(df$Species, levels = c("AA", "CQ", "CP"))

new_location <- c("Field - Guadeloupe", 
                  "Laboratory - Slab TC (Wolbachia -)",
                  "Laboratory - Lavar",
                  "Field - Bosc", 
                  "Field - Camping Europe")
df$Strain <- factor(df$Strain, levels = new_location)

levels(df$Strain)

## [1] "Field - Guadeloupe"                 "Laboratory - Slab TC (Wolbachia -)"
## [3] "Laboratory - Lavar"                 "Field - Bosc"                      
## [5] "Field - Camping Europe"

# convert df into data.table and set min and max for y-axis
dat <- data.table::data.table(df)
dat[, y_min := value*0.1, by = variable]
dat[, y_max := value*1.1 , by = variable]
dat[dat$variable=="Observed", "y_min"] <- 10
dat[dat$variable=="Chao1", "y_min"] <- 10
dat[dat$variable=="Observed", "y_max"] <- 55
dat[dat$variable=="Chao1", "y_max"] <- 55
#dat[dat$variable=="Shannon", "y_max"] <- 3

# plot
p <- ggplot(dat,aes(Species,value, color=Strain)) +
  facet_wrap(~ variable, drop=T, scale="free_y")+
  geom_boxplot(outlier.colour = NA, alpha=0.8, position = position_dodge2(width=1, preserve="single"))+
  geom_point(size=1, aes(shape=dat$Strain), position=position_dodge(width=0.7, preserve='total'))+
  scale_shape_manual("Strain", values = c(16,17,17,16,16),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), " -)")),
                                "Laboratory - Lavar",
                                "Field - Bosc", 
                                "Field - Camping Europe"))+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#5B6BF7", "#00B0F6", "#A3A500", "#F8766D"),
                     labels = c("Field - Guadeloupe", 
                                expression(paste("Laboratory - Slab TC (", italic("Wolbachia"), " -)")),
                                "Laboratory - Lavar",
                                "Field - Bosc", 
                                "Field - Camping Europe"))+
  labs(x="Species", y = "Alpha Diversity Measure")+
  theme(legend.text.align = 0)
p

# final plot with more space on y-axis
p_bis <- p +
  geom_blank(aes(y = y_min)) +
  geom_blank(aes(y = y_max))
p_bis

Whole with rarefaction

# Rarefaction with sample.size=100
ps.whole2 <- ps.whole %>% rarefy_even_depth(rngseed=1, sample.size = 100)

## `set.seed(1)` was used to initialize repeatable random subsampling.

## Please record this for your records so others can reproduce.

## Try `set.seed(1); .Random.seed` for the full vector

## ...

## 1 samples removedbecause they contained fewer reads than `sample.size`.

## Up to first five removed samples are:

## NP34

## ...

## 4OTUs were removed because they are no longer 
## present in any sample after random subsampling

## ...

ps.whole2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 63 taxa and 66 samples ]
## sample_data() Sample Data:       [ 66 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 63 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 63 reference sequences ]

ps.whole2@sam_data %>% data.frame() -> test2
compute_read_counts(ps.whole2)

## [1] 6600

p2 <- alpha_div(ps.whole2, measures)
p2

## Rarefaction with sample.size=1000
ps.whole3 <- ps.whole %>% rarefy_even_depth(rngseed=2, sample.size = 1000)

## `set.seed(2)` was used to initialize repeatable random subsampling.

## Please record this for your records so others can reproduce.

## Try `set.seed(2); .Random.seed` for the full vector

## ...

## 4 samples removedbecause they contained fewer reads than `sample.size`.

## Up to first five removed samples are:

## NP29NP30NP34NP36

## ...

ps.whole3

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 67 taxa and 63 samples ]
## sample_data() Sample Data:       [ 63 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 67 taxa by 7 taxonomic ranks ]
## refseq()      DNAStringSet:      [ 67 reference sequences ]

ps.whole3@sam_data %>% data.frame() -> test3
compute_read_counts(ps.whole3)

## [1] 63000

p3 <- alpha_div(ps.whole3, measures)
p3

Ovaries

# select ovary samples
ps.ovary <- subset_samples(ps.filter, Organ=="Ovary")
sample_data(ps.ovary)$Read_depth <- sample_sums(ps.ovary)

# change factors and order in metadata
metadata.ps <- data.frame(sample_data(ps.ovary))
levels(metadata.ps$Strain)

## [1] "Field - Bosc"           "Field - Camping Europe" "Field - Guadeloupe"    
## [4] "Laboratory - Lavar"

levels(metadata.ps$Strain) <- c(levels(metadata.ps$Strain), "Laboratory - Lavar")
metadata.ps$Strain[metadata.ps$Strain=="Lavar (labo)"] <- "Laboratory - Lavar"
levels(metadata.ps$Strain)

## [1] "Field - Bosc"           "Field - Camping Europe" "Field - Guadeloupe"    
## [4] "Laboratory - Lavar"

metadata.ps$Strain <- droplevels(metadata.ps$Strain)
levels(metadata.ps$Strain)

## [1] "Field - Bosc"           "Field - Camping Europe" "Field - Guadeloupe"    
## [4] "Laboratory - Lavar"

# update metadata in phyloseq object
sample_data(ps.ovary) <- metadata.ps

# estimate alpha diversity
p1 <- plot_richness(ps.ovary, 
                      x="Sample", 
                      color="Strain", 
                      measures=measures, 
                      nrow = 1)

# extract data from richness plot to custom it
df <- p1$data

# changer levels and order of species and location for plot
levels(df$Species) <- c("AA", "CP", "CQ")

df$Species <- factor(df$Species, levels = c("AA", "CQ", "CP"))

new_location <- c("Field - Guadeloupe", 
                  "Laboratory - Slab TC (Wolbachia -)",
                  "Field - Bosc", 
                  "Field - Camping Europe", 
                  "Laboratory - Lavar")
df$Strain <- factor(df$Strain, levels = new_location)

# final plot
p_ovary <- ggplot(df,aes(Species,value,colour=Strain)) +
  facet_wrap(~ variable, drop=T,scale="free")+
  scale_color_manual("Strain", 
                     values=c("#00BF7D", "#A3A500", "#F8766D", "#00B0F6"),
                     labels = c("Field - Guadeloupe", 
                                "Field - Bosc", 
                                "Field - Camping Europe", 
                                "Laboratory - Lavar"))+
  geom_boxplot(outlier.colour = NA, alpha=0.8, position = position_dodge2(width=1, preserve="single")) +
  geom_point(size=2, position=position_dodge(width=0.7, preserve='total'))+
  labs(x="Species", y = "Alpha Diversity Measure")+
    theme(legend.text.align = 0)

p_ovary

Mean of estimated alpha diversity

df_score <- data.frame(p$data)

x <- c("Culex pipiens - Field", "Culex pipiens - Bosc", "Culex pipiens - Camping Europe", "Culex pipiens - Lavar (lab)", 
       "Culex quinquefasciatus - Guadeloupe", "Culex quinquefasciatus - Slab TC",
       "Aedes aegypti - Guadeloupe")
y <- c("Observed", "Chao1", "Shannon")

df <- data.frame(matrix(ncol=3, nrow=7))
rownames(df) <- x
colnames(df) <- y

# Aedes aegypti (from Guadeloupe)
df_score[df_score$Species=="AA" & df_score$variable=="Observed", "value"] %>% mean() -> df[7,1]
df_score[df_score$Species=="AA" & df_score$variable=="Chao1", "value"] %>% mean() -> df[7,2]
df_score[df_score$Species=="AA" & df_score$variable=="Shannon", "value"] %>% mean() -> df[7,3]

# Culex pipiens from field
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Field=="Field", "value"] %>% mean() -> df[1,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Field=="Field", "value"] %>% mean() -> df[1,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Field=="Field", "value"] %>% mean() -> df[1,3]

# Culex pipiens from labo
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Laboratory - Lavar", "value"] %>% mean() -> df[4,3]

# Culex pipiens from Bosc
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Field - Bosc", "value"] %>% mean() -> df[2,3]

# Culex pipiens from Camping Europe
df_score[df_score$Species=="CP" & df_score$variable=="Observed" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,1]
df_score[df_score$Species=="CP" & df_score$variable=="Chao1" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,2]
df_score[df_score$Species=="CP" & df_score$variable=="Shannon" & df_score$Strain=="Field - Camping Europe", "value"] %>% mean() -> df[3,3]

# Culex quinquefasciatus from field (Guadeloupe)
df_score[df_score$Species=="CQ" & df_score$variable=="Observed" & df_score$Field=="Field", "value"] %>% mean() -> df[5,1]
df_score[df_score$Species=="CQ" & df_score$variable=="Chao1" & df_score$Field=="Field", "value"] %>% mean() -> df[5,2]
df_score[df_score$Species=="CQ" & df_score$variable=="Shannon" & df_score$Field=="Field", "value"] %>% mean() -> df[5,3]

# Culex quinquefasciatus from labo (Wolbachia-)
df_score[df_score$Species=="CQ" & df_score$variable=="Observed" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,1]
df_score[df_score$Species=="CQ" & df_score$variable=="Chao1" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,2]
df_score[df_score$Species=="CQ" & df_score$variable=="Shannon" & df_score$Field=="Lab ", "value"] %>% mean() -> df[6,3]

df %>%
  kbl(booktable=TRUE) %>%
  kable_paper("hover", full_width = F)

	Observed	Chao1	Shannon
Culex pipiens - Field	20.84211	22.17982	0.1971403
Culex pipiens - Bosc	20.61538	20.94231	0.1790446
Culex pipiens - Camping Europe	21.33333	24.86111	0.2363475
Culex pipiens - Lavar (lab)	27.90909	30.03485	1.0724543
Culex quinquefasciatus - Guadeloupe	17.25000	24.66667	1.5301268
Culex quinquefasciatus - Slab TC	28.46154	31.08269	1.8549403
Aedes aegypti - Guadeloupe	32.22222	33.51852	1.0625983

Save plots

tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.tiff", units="in", width=8, height=5, res=300)
p
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.png", units="in", width=8, height=5, res=300)
p
dev.off()

## quartz_off_screen 
##                 2

pdf ( "../../../../output/1_MED/1E/1Eb_MED_alpha_diversity.pdf")
p_bis
dev.off()

## quartz_off_screen 
##                 2

tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_100.tiff", units="in", width=8, height=5, res=300)
p2
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_100.png", units="in", width=8, height=5, res=300)
p2
dev.off()

## quartz_off_screen 
##                 2

tiff("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_1000.tiff", units="in", width=8, height=5, res=300)
p3
dev.off()

## quartz_off_screen 
##                 2

png("../../../../output/1_MED/1E/1Eb_MED_alpha_diversity_rare_1000.png", units="in", width=8, height=5, res=300)
p3
dev.off()

## quartz_off_screen 
##                 2

Statistics

Prepare subdatasets for stats

# Whole
df <- p$data
df_culex <- df[df$Species!="AA",]
df_pipiens <- df[df$Species=="CP",]
df_quinque <- df[df$Species=="CQ",]

# Whole rarefied (100)
df2 <- p2$data
df2_culex <- df2[df2$Species!="AA",]
df2_pipiens <- df2[df2$Species=="CP",]
df2_quinque <- df2[df2$Species=="CQ",]

# Whole rarefied (1000)
df3 <- p3$data
df3_culex <- df3[df3$Species!="AA",]
df3_pipiens <- df3[df3$Species=="CP",]
df3_quinque <- df3[df3$Species=="CQ",]

# Wolbachia- vs Wolbachia+
df_culex_wolbachia <- df_culex[df_culex$Species=="CQ" & df_culex$Strain!="Laboratory - Lavar",]

# Wolbachia- vs Wolbachia+ (rarefied 100)
df2_culex_wolbachia <- df2_culex[df_culex$Species=="CQ" & df2_culex$Strain!="Laboratory - Lavar",]

# Wolbachia- vs Wolbachia+ (rarefied 1000)
df3_culex_wolbachia <- df3_culex[df_culex$Species=="CQ" & df3_culex$Strain!="Laboratory - Lavar",]

Influence of depth sequencing

Effect of Species

Observed

# names for final array
names <- c("Group", "Read_depth", "Group:Read_depth", "Residuals")
names2 <- c("Group", "Read_depth", "Residuals")

## Observed
df_observed <- df[df$variable=="Observed"  & df$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
wb <- createWorkbook()
addWorksheet(wb, "Alpha diversity (depth seq)")

writeData(wb, "Alpha diversity (depth seq)", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 6, rowNames = TRUE)

## Observed (rarefied 100)
df_observed_100 <- df2[df$variable=="Observed"  & df2$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 16, rowNames = TRUE)

## Observed (rarefied, 1000)
df_observed_1000 <- df3[df3$variable=="Observed"  & df3$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 26, rowNames = TRUE)

Chao1

# Chao1
df_chao <- df[df$variable=="Chao1"  & df$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 6, rowNames = TRUE)

# Chao1 (rarefied, 100)
df_chao_100 <- df2[df2$variable=="Chao1"  & df2$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 16, rowNames = TRUE)

## Chao (rarefied, 1000)
df_chao_1000 <- df3[df3$variable=="Chao1"  & df3$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 1000)", startCol = 9, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 26, rowNames = TRUE)
#saveWorkbook(wb, "Supplementary_Table_2_review_sheet1_30_11_21.xlsx", overwrite = TRUE)

Shannon

## Shannon
df_shannon <- df[df$variable=="Shannon"  & df$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 3)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 5)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 6, rowNames = TRUE)

## Shannon (rarefied, 100)
df_shannon_100 <- df2[df2$variable=="Shannon"  & df2$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 13)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 15)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 16, rowNames = TRUE)

## Shannon (rarefied, 1000)
df_shannon_1000 <- df3[df3$variable=="Shannon"  & df3$Field!="Lab ",]

#### Construct linear model
model  <- lm(value~ Species*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 23)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 25)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 26, rowNames = TRUE)

saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of tetracycline

Observed

## Observed
df_observed<- df_quinque[df_quinque$variable=="Observed",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 35)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 40, rowNames = TRUE)

## Observed (rarefied 100)
df_observed_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Observed",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 50, rowNames = TRUE)


## Observed (rarefied, 1000)
df_observed_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Observed",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 60, rowNames = TRUE)

Chao1

## Chao1
df_chao <- df_quinque[df_quinque$variable=="Chao1",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 40, rowNames = TRUE)

## Chao1 (rarefied 100)
df_chao_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Chao1",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 50, rowNames = TRUE)


## Chao1 (rarefied, 1000)
df_chao_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Chao1",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
#ggpubr::ggqqplot(residuals(model))
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO INDEX (rarefied, 1000)", startCol = 9, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 60, rowNames = TRUE)

Shannon

## Shannon
df_shannon <- df_quinque[df_quinque$variable=="Shannon",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
#ggpubr::ggqqplot(residuals(model))
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 37)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 39)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 40, rowNames = TRUE)

## Shannon (rarefied 100)
df_shannon_100 <- df2_culex_wolbachia[df2_culex_wolbachia$variable=="Shannon",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 47)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 49)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 50, rowNames = TRUE)

## Shannon (rarefied, 1000)
df_shannon_1000 <- df3_culex_wolbachia[df3_culex_wolbachia$variable=="Shannon",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 57)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 59)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 60, rowNames = TRUE)

saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)

Effect of laboratory

Observed

## Observed
df_observed <- df_pipiens[df_pipiens$variable=="Observed",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "Effect of lab and location in the Culex pipiens samples", startCol = 1, startRow = 68)
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX", startCol = 1, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 72, rowNames = TRUE)

## Observed (rarefied 100)
df_observed_100 <- df2_pipiens[df2_pipiens$variable=="Observed",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 100)", startCol = 1, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 82, rowNames = TRUE)

## Observed (rarefied, 1000)
df_observed_1000 <- df3_pipiens[df3_pipiens$variable=="Observed", ]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_observed_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "OBSERVED INDEX (rarefied, 1000)", startCol = 1, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 1, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 1, startRow = 92, rowNames = TRUE)

Chao1

## Chao1
df_chao <- df_pipiens[df_pipiens$variable=="Chao1",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX", startCol = 9, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 72, rowNames = TRUE)

## Chao1 (rarefied 100)
df_chao_100 <- df2_pipiens[df2_pipiens$variable=="Chao1",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

writeData(wb, "Alpha diversity (depth seq)", "CHAO1 INDEX (rarefied, 100)", startCol = 9, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 82, rowNames = TRUE)

## Chao1 (rarefied, 1000)
df_chao_1000 <- df3_pipiens[df3_pipiens$variable=="Chao1", ]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_chao_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "CHAO INDEX (rarefied, 1000)", startCol = 9, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 9, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 9, startRow = 92, rowNames = TRUE)

Shannon

## Shannon
df_shannon <- df_pipiens[df_pipiens$variable=="Shannon",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX", startCol = 17, startRow = 70)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 71)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 72, rowNames = TRUE)


## Shannon (rarefied 100)
df_shannon_100 <- df2_pipiens[df2_pipiens$variable=="Shannon",]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon_100)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 100)", startCol = 17, startRow = 79)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 81)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 82, rowNames = TRUE)

## Shannon (rarefied, 1000)
df_shannon_1000 <- df3_pipiens[df3_pipiens$variable=="Shannon", ]

#### Construct linear model
model  <- lm(value~ Strain*Read_depth, data = df_shannon_1000)
anova(model)-> res_anova
res_anova <- res_anova %>% fill_significance("Pr(>F)")

#### Create residuals plot
par(mfrow=c(2,2))
plot(model)

#### Excel
writeData(wb, "Alpha diversity (depth seq)", "SHANNON INDEX (rarefied, 1000)", startCol = 17, startRow = 89)
writeData(wb, "Alpha diversity (depth seq)", "ANOVA:", startCol = 17, startRow = 91)
writeData(wb, "Alpha diversity (depth seq)", res_anova, startCol = 17, startRow = 92, rowNames = TRUE)

saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)

Analysis of alpha diversity

Effect of species

Observed

names <- c("Species", "Residuals")

## Observed
df_observed <- df[df$variable=="Observed"  & df$Field!="Lab ",]
cat("OBSERVED\n\n")

## OBSERVED

res <- make_anova_tukey(df_observed, df_observed$Species, df_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## var        2  974.89  487.44  10.543 0.0003622 ***
## Residuals 29 1340.83   46.24                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## CQ - AA == 0  -14.972      4.086  -3.664  0.00265 ** 
## CP - AA == 0  -11.380      2.751  -4.136  < 0.001 ***
## CP - CQ == 0    3.592      3.741   0.960  0.60166    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
addWorksheet(wb, "Alpha diversity")

writeData(wb, "Alpha diversity", "Effect of species in the field samples", startCol = 1, startRow = 1)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 12, rowNames = TRUE)

Chao1

## Chao1
df_chao <- df[df$variable=="Chao1"  & df$Field!="Lab ",]
cat("CHAO1\n\n")

## CHAO1

res <- make_anova_tukey(df_chao, df_chao$Species, df_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df  Sum Sq Mean Sq F value   Pr(>F)   
## var        2  789.86  394.93  8.5122 0.001234 **
## Residuals 29 1345.47   46.40                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## CQ - AA == 0   -8.852      4.093  -2.163   0.0923 .  
## CP - AA == 0  -11.339      2.756  -4.114   <0.001 ***
## CP - CQ == 0   -2.487      3.747  -0.664   0.7824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 12, rowNames = TRUE)

Shannon

## Shannon
df_shannon <- df[df$variable=="Shannon"  & df$Field!="Lab ",]
cat("SHANNON\n\n")

## SHANNON

res <- make_anova_tukey(df_shannon, df_shannon$Species, df_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## var        2 8.4685  4.2342  46.442 9.09e-10 ***
## Residuals 29 2.6440  0.0912                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##              Estimate Std. Error t value Pr(>|t|)    
## CQ - AA == 0   0.4675     0.1814   2.577   0.0384 *  
## CP - AA == 0  -0.8655     0.1222  -7.083   <0.001 ***
## CP - CQ == 0  -1.3330     0.1661  -8.025   <0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 3)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 5)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 6, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 11)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 12, rowNames = TRUE)

Effect of tetracycline

names <- c("Strain", "Residuals")

## Observed
culex_observed <- df_quinque[df_quinque$variable=="Observed",]
cat("OBSERVED\n\n")

## OBSERVED

res <- make_anova_tukey(culex_observed, culex_observed$Strain, culex_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## var        1 384.49  384.49  21.521 0.0003211 ***
## Residuals 15 267.98   17.87                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                              Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0   11.212
##                                                              Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0      2.417   4.639
##                                                              Pr(>|t|)    
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0 0.000321 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "Effect of tetracycline in the Culex quinquefasciatus samples", startCol = 1, startRow = 20)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 31, rowNames = TRUE)

Chao1

## Chao1
culex_chao <- df_quinque[df_quinque$variable=="Chao1",]
cat("CHAO1\n\n")

## CHAO1

res <- make_anova_tukey(culex_chao, culex_chao$Strain, culex_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value  Pr(>F)  
## var        1 125.92 125.918  5.3198 0.03577 *
## Residuals 15 355.04  23.669                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                              Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0    6.416
##                                                              Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0      2.782   2.306
##                                                              Pr(>|t|)  
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0   0.0358 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 31, rowNames = TRUE)

Shannon

## Shannon
culex_shannon <- df_quinque[df_quinque$variable=="Shannon",]
cat("SHANNON\n\n")

## SHANNON

res <- make_anova_tukey(culex_shannon, culex_shannon$Strain, culex_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df  Sum Sq Mean Sq F value Pr(>F)
## var        1 0.32272 0.32272  1.7577 0.2047
## Residuals 15 2.75404 0.18360               
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                              Estimate
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0   0.3248
##                                                              Std. Error t value
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0     0.2450   1.326
##                                                              Pr(>|t|)
## Laboratory - Slab TC (Wolbachia -) - Field - Guadeloupe == 0    0.205
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 22)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 24)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 25, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 30)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 31, rowNames = TRUE)

Effect of lab in Culex pipiens samples

Observed

## Observed
pipiens_observed <- df_pipiens[df_pipiens$variable=="Observed",]
cat("OBSERVED\n\n")

## OBSERVED

res <- make_anova_tukey(pipiens_observed, pipiens_observed$Strain, pipiens_observed$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## var        2 511.28 255.642  13.195 4.45e-05 ***
## Residuals 38 736.23  19.374                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                  Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0            -7.2937     1.5398  -4.737
## Field - Camping Europe - Laboratory - Lavar == 0  -6.5758     2.0272  -3.244
## Field - Camping Europe - Field - Bosc == 0         0.7179     2.1724   0.330
##                                                  Pr(>|t|)    
## Field - Bosc - Laboratory - Lavar == 0            < 0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0  0.00682 ** 
## Field - Camping Europe - Field - Bosc == 0        0.94077    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "Effect of lab (and location) in the Culex pipiens samples", startCol = 1, startRow = 37)
writeData(wb, "Alpha diversity", "OBSERVED INDEX", startCol = 1, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 1, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 1, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 1, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 1, startRow = 47, rowNames = TRUE)

Chao1

## Chao1
pipiens_chao <- df_pipiens[df_pipiens$variable=="Chao1",]
cat("CHAO1\n\n")

## CHAO1

res <- make_anova_tukey(pipiens_chao, pipiens_chao$Strain, pipiens_chao$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## var        2  692.1  346.05  10.453 0.0002415 ***
## Residuals 38 1258.0   33.11                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                  Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0             -9.093      2.013  -4.517
## Field - Camping Europe - Laboratory - Lavar == 0   -5.174      2.650  -1.952
## Field - Camping Europe - Field - Bosc == 0          3.919      2.840   1.380
##                                                  Pr(>|t|)    
## Field - Bosc - Laboratory - Lavar == 0             <0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0    0.135    
## Field - Camping Europe - Field - Bosc == 0          0.357    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "CHAO1 INDEX", startCol = 9, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 9, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 9, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 9, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 9, startRow = 47, rowNames = TRUE)

Shannon

## Shannon
pipiens_shannon <- df_pipiens[df_pipiens$variable=="Shannon",]
cat("SHANNON\n\n")

## SHANNON

res <- make_anova_tukey(pipiens_shannon, pipiens_shannon$Strain, pipiens_shannon$value)
res_tukey <- glht.table(res[[2]]) %>% as.data.frame()
res_tukey <- fill_significance(res_tukey, "Pr(>|t|)")
res_anova <- res[[1]] %>% as.data.frame()
res_anova <- fill_significance(res_anova,"Pr(>F)")
rownames(res_anova) <- names
res

## $anova
## Analysis of Variance Table
## 
## Response: value
##           Df Sum Sq Mean Sq F value    Pr(>F)    
## var        2 7.8247  3.9124  18.859 2.047e-06 ***
## Residuals 38 7.8832  0.2075                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $tukey
## 
##   Simultaneous Tests for General Linear Hypotheses
## 
## Multiple Comparisons of Means: Tukey Contrasts
## 
## 
## Fit: lm(formula = value ~ var, data = df)
## 
## Linear Hypotheses:
##                                                  Estimate Std. Error t value
## Field - Bosc - Laboratory - Lavar == 0            -0.8934     0.1593  -5.607
## Field - Camping Europe - Laboratory - Lavar == 0  -0.8361     0.2098  -3.986
## Field - Camping Europe - Field - Bosc == 0         0.0573     0.2248   0.255
##                                                  Pr(>|t|)    
## Field - Bosc - Laboratory - Lavar == 0             <0.001 ***
## Field - Camping Europe - Laboratory - Lavar == 0   <0.001 ***
## Field - Camping Europe - Field - Bosc == 0          0.964    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)

### Excel
writeData(wb, "Alpha diversity", "SHANNON INDEX", startCol = 17, startRow = 39)
writeData(wb, "Alpha diversity", "ANOVA:", startCol = 17, startRow = 41)
writeData(wb, "Alpha diversity", res_anova, startCol = 17, startRow = 42, rowNames = TRUE)
writeData(wb, "Alpha diversity", "Tukey:", startCol = 17, startRow = 46)
writeData(wb, "Alpha diversity", res_tukey, startCol = 17, startRow = 47, rowNames = TRUE)

saveWorkbook(wb, "../../../../output/1_MED/1E/Supplementary_Table_2.xlsx", overwrite = TRUE)

Session info

sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.3.1 openxlsx_4.2.3   multcomp_1.4-14  TH.data_1.0-10  
##  [5] MASS_7.3-53      survival_3.2-7   mvtnorm_1.1-1    forcats_0.5.0   
##  [9] stringr_1.4.0    dplyr_1.0.2      purrr_0.3.4      readr_1.4.0     
## [13] tidyr_1.1.2      tibble_3.0.4     tidyverse_1.3.0  ggplot2_3.3.2   
## [17] phyloseq_1.30.0 
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-149        fs_1.5.0            lubridate_1.7.9    
##  [4] webshot_0.5.2       httr_1.4.2          tools_3.6.3        
##  [7] backports_1.1.10    bslib_0.2.4         R6_2.4.1           
## [10] vegan_2.5-6         DBI_1.1.0           BiocGenerics_0.32.0
## [13] mgcv_1.8-33         colorspace_2.0-0    permute_0.9-5      
## [16] ade4_1.7-15         withr_2.3.0         tidyselect_1.1.0   
## [19] compiler_3.6.3      cli_2.1.0           rvest_0.3.6        
## [22] Biobase_2.46.0      xml2_1.3.2          sandwich_3.0-0     
## [25] labeling_0.4.2      bookdown_0.22       sass_0.3.1         
## [28] scales_1.1.1        digest_0.6.26       rmarkdown_2.7      
## [31] XVector_0.26.0      pkgconfig_2.0.3     htmltools_0.5.1.1  
## [34] highr_0.8           dbplyr_1.4.4        rlang_0.4.10       
## [37] readxl_1.3.1        rstudioapi_0.11     farver_2.0.3       
## [40] jquerylib_0.1.3     generics_0.0.2      zoo_1.8-8          
## [43] jsonlite_1.7.1      zip_2.1.1           magrittr_1.5       
## [46] biomformat_1.14.0   Matrix_1.2-18       fansi_0.4.1        
## [49] Rcpp_1.0.5          munsell_0.5.0       S4Vectors_0.24.4   
## [52] Rhdf5lib_1.8.0      ape_5.4-1           lifecycle_0.2.0    
## [55] stringi_1.5.3       yaml_2.2.1          zlibbioc_1.32.0    
## [58] rhdf5_2.30.1        plyr_1.8.6          grid_3.6.3         
## [61] blob_1.2.1          parallel_3.6.3      crayon_1.3.4       
## [64] lattice_0.20-41     Biostrings_2.54.0   haven_2.3.1        
## [67] splines_3.6.3       multtest_2.42.0     hms_0.5.3          
## [70] ps_1.4.0            knitr_1.30          pillar_1.4.6       
## [73] igraph_1.2.6        reshape2_1.4.4      codetools_0.2-16   
## [76] stats4_3.6.3        reprex_0.3.0        glue_1.4.2         
## [79] evaluate_0.14       data.table_1.13.2   modelr_0.1.8       
## [82] vctrs_0.3.4         rmdformats_1.0.2    foreach_1.5.1      
## [85] cellranger_1.1.0    gtable_0.3.0        assertthat_0.2.1   
## [88] xfun_0.22           broom_0.7.2         viridisLite_0.3.0  
## [91] iterators_1.0.13    IRanges_2.20.2      cluster_2.1.0      
## [94] ellipsis_0.3.1